data = data3d('../data/head.npz', 'short')
cut = float_param(-1.8, -1.8, 1.8, 'C')

bone_alpha = float_param(0.1, 0, 0.5, 'B')
soft_alpha = float_param(0.1, 0, 0.5, 'S')
water_alpha = float_param(0.1, 0, 0.5, 'W')
fat_alpha = float_param(0.1, 0, 0.5, 'F')

density = cubic_query_3d_cut(data, X, D, cut) * 32768
rgba = float4(0)

if density > 200:
    rgba = float4(1, 1, 1, bone_alpha)
elif density > 50:
    rgba = float4(0.5, 0.3, 0.3, soft_alpha)
elif density > -50:
    rgba = float4(0.1, 0.1, 0.5, water_alpha)
elif density > -100:
    rgba = float4(0.5, 0.5, 0.1, fat_alpha)
else:
    pass

if rgba.w == 0: return 0     

# light direction (you can make this tweakable too)
L = normalize(float3(-2, -2, -2) - X)

# compute phong lighting
N = -normalize(cubic_gradient_3d_cut(data, X, D, cut))
color = phong(L, N, -D, rgba.xyz, 0.8, 50, 0.2)

return float4(color.x, color.y, color.z, rgba.w)

